Aisha: Não tem o que fazer.
Aisha: É para usar o filtro automático mas explicar bem o que ele faz.
4390- Taxa de juros (Selic acumulada no mês). Qual a diferença entre essas duas?Aisha: Usar a SELIC do Bacen é melhor (depois tem que voltar aqui e escrever quais as diferenças entre essas séries e porque a do Bacen é melhor), porém tem que lembrar de passar para valores anuais (por enquanto estava mensal). Update: Fui procurar e entendi o que é a Selic Overnight (que tem no IPEA), mas ainda não entendi o que é a Selic acumulada no mês do Banco Central. Update: A Selic Acumulada no site do Bacen e a Selic Overnight do site do IPEA tem os mesmos valores. E a Selic Acumulada é feita usando a Selic diária (série 11) com a fórmula de juros compostos. Agora resta saber por que a série 11 tem formato tão de escadinha.
Aisha: Olhar o que o Mumtaz usou no VAR dele. Basicamente procurar por produto, inflação, câmbio (talvez)…
Aisha: De acordo com o Guilherme, se usar câmbio o melhor é fazer a partir do fim da âncora cambial. Fui olhar e isso implica começar a série em Janeiro/99, quando foi adotado o regime de metas de inflação (http://www.rep.org.br/pdf/87-1.pdf). Por enquanto eu vou rodar desde 94 e depois eu adapto, se for o caso. Fui ver e a série do IBC-Br começa em 2003, daí fode tudo.
Aisha: Dá para usar isso? https://cran.r-project.org/web/packages/GetTDData/vignettes/gtdd-vignette_GetTDData.html . Update: O Portela disse o seguinte " Oi Aisha, Você poderia usar a taxa do swap DI-Pré de 3 meses. Dá pra baixar do banco central ou então pelo BETS. Um abraço“. No banco central, a série é 7818- Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período). Mas eu não entendi que trem é esse… como que isso se relaciona com política monetária? E qual seria a justificativa para usar isso e não a taxa de juros da economia para colocar no VAR do Mumtaz?
Aisha: Dar uma olhada no material de macro, uma delas era taxa de crescimento (diferença dos logs).
Aisha: Sei lá. O Mumtaz usa growth of the nominal effective exchange rate, mas aí não sei se ele faz a diferença dos logs e por isso chama assim ou se ele calcula a taxa de crescimento e depois faz a diferença dos logs. Eu posso dar uma olhada no código deles depois.
library(ggplot2, quietly = TRUE)
library(forecast, quietly = TRUE)
library(BETS, quietly = TRUE)
library(seasonal, quietly = TRUE)
library(seasonalview, quietly = TRUE)
library(lubridate, quietly = TRUE)
library(zoo, quietly = TRUE)
library(stargazer, quietly = TRUE)
library(gridExtra, quietly = TRUE)
library(reshape2, quietly = TRUE)
library(ggfortify, quietly = TRUE)
library(RColorBrewer, quietly = TRUE)
library(scales, quietly = TRUE)
library(quantmod, quietly = TRUE)
library(PerformanceAnalytics, quietly = TRUE)
library(strucchange, quietly = TRUE)
Aisha: Caso o BETS não esteja instalado, retirar o eval=FALSE do chunk abaixo e rodar ele.
library(devtools)
devtools::install_github("pedrocostaferreira/BETS", force= TRUE)
Os dados podem ser pegos de https://www3.bcb.gov.br/sgspub/localizarseries/localizarSeries.do?method=prepararTelaLocalizarSeries:
Renda do trabalho: Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (7620)
Renda do capital: Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (7621)
Ambas séries são mensais e começam em 31/01/1992. O último dado disponível em 05/01/2018 era de novembro de 2017, totalizando 311 meses. O pacote BETS consegue fazer download das séries apenas usando o código delas. Cortando em 1994 (por causa do plano real) ficam 281 observações (vou deixar para ver o trem da âncora cambial depois).
Update: Por enquanto estou incluindo IBC-Br o que me deixa com 177 observações apenas…
trabalho <- BETS.get("7620")
capital <- BETS.get("7620")
autoplot(trabalho)
autoplot(capital)
trabalho <- BETS.get("7620", from = "1994-07-01", to = "2017-11-01")
autoplot(trabalho)
capital <- BETS.get("7621", from = "1994-07-01", to = "2017-11-01")
autoplot(capital)
capital_trabalho <- capital/trabalho
autoplot(capital/trabalho)
length(capital_trabalho)
## [1] 281
Inicialmente eu fiz um primeiro VAR apenas com Selic e razão capital/trabalho. Havia a dúvida se era melhor usar a Selic do site do BACEN ou a Selic Overnight do site do IPEA. De acordo com o Guilherme, o melhor é o do Bacen pois é a efetiva (não entendi bem essa história…), mas lá em Bayesiana ele não me disse isso “porque estava ocupado me ensinando o que era Selic”.
Do site do Banco Central, tem essa definição: * “Define-se Taxa Selic como a taxa média ajustada dos financiamentos diários apurados no Sistema Especial de Liquidação e de Custódia (Selic) para títulos federais. Para fins de cálculo da taxa, são considerados os financiamentos diários relativos às operações registradas e liquidadas no próprio Selic e em sistemas operados por câmaras ou prestadores de serviços de compensação e de liquidação” (Link).
Aqui tem uma explicação do que é a Selic Overnight:
Mas isso ainda não explica o que é a série 4390 - Taxa de juros - Selic acumulada no mês que tem no sistema do Banco Central. Olhando na caixinha de explicações metodológicas, diz que ela é derivada da série 11 - Taxa de Juros - Selic. Essa taxa 11 está em %a.d.. Então vou pegar uns valores para ver o que acontece:
selic_11 <- BETS.get("11", from = "2000-01-01", to = "2017-12-31")
#head(selic_11)
selic_4390 <- BETS.get("4390", from = "2000-01-01", to = "2017-12-31")
#head(selic_4390)
df1 <- as.data.frame(selic_11)
#head(df1)
p1 <- ggplot(df1, aes(date, value)) +
geom_line() +
scale_x_date() +
xlab("Data") +
ylab("Selic (%a.d.)") +
ggtitle("11 - Taxa de Juros - Selic")
df2 <- data.frame(time(selic_4390), as.numeric(selic_4390))
#head(df2)
names(df2) <- c("date", "value")
p2 <- ggplot(df2, aes(date, value)) +
geom_line() +
xlab("Data") +
ylab("Selic (%a.m.)") +
ggtitle("4390 - Taxa de Juros - Selic acumulada no mês")
grid.arrange(p1, p2, ncol=1, nrow = 2)
## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
Depois disso, eu peguei ambas séries de janeiro de 2017 a dezembro de 2017 e na série diária eu somei os valores de cada mês, para obter 12 valores e comparei com os 12 valores da série 4390. Os valores são bem similares. Depois eu tentei fazer um puxadinho do que é a fórmula nos dados básicos da série 4390, onde diz CONVPER((ACMVALORES(((SERIE(11)/100)+1),"mensal","multiplicacao")-1)*100,"mensal","ultimovalor"). O que eu fiz no Excel foi pegar cada valor, dividir por \(100\) e somar \(1\). Depois disso, multipliquei os valores dentro de um mesmo mês, subtraí \(1\) e multipliquei por \(100\).
Comparação da série 11 com a série 4390
Resumindo, o que tem na série 4390 parece ser a fórmula de juros compostos usando os valores da série 11, que por sua vez tem valores do tipo “escadinha”.
Agora, vou olhar para a Selic Overnight, do site do IPEA.
selic_overnight <- read.csv2("C:\\Users\\Aishameriane\\OneDrive\\Documentos\\Mestrado Economia\\Teoria Macroeconômica II - 2017-2\\Artigo\\Aplicação\\selic_overnight.csv", sep = ";", dec = ",", header = TRUE)
selic_overnight <- selic_overnight[,-3]
colnames(selic_overnight) <- c("Data", "Selic")
# Arrumando as datas
selic_overnight[,1]<-paste(selic_overnight[,1], ".01", sep="")
selic_overnight[,1]<-ymd(selic_overnight[,1])
inicio <- which(selic_overnight[,1] == "1994-07-01")
fim <- which(selic_overnight[,1] == "2017-11-01")
selic_overnight<-selic_overnight[inicio:fim,]
p3 <- ggplot(selic_overnight, aes(Data, Selic)) +
geom_line() +
xlab("Data") +
ylab("Selic (%a.m.)") +
ggtitle("Taxa de Juros - Overnight Selic")
grid.arrange(p1, p2, p3, ncol=1, nrow = 3)
Eu estou desconfiada que a Overnight e a série 4390 são a mesma coisa, então vou dropar uns valores da overnight e ver o que acontece.
inicio <- which(selic_overnight[,1] == "2000-01-01")
fim <- which(selic_overnight[,1] == "2017-11-01")
selic_overnight_menor<-selic_overnight[inicio:fim,]
p4 <- ggplot(selic_overnight_menor, aes(Data, Selic)) +
geom_line() +
xlab("Data") +
ylab("Selic (%a.m.)") +
ggtitle("Taxa de Juros - Overnight Selic")
grid.arrange(p1, p2, p4, ncol=1, nrow = 3)
Os gráficos são bem similares… Mas e os valores? (obs: eu olhei para os valores da cauda também e usei janelas maiores de dados… só os dois últimos não são iguais na segunda casa decimal.)
cbind(head(round(selic_overnight_menor[,2],2), 10), head(selic_4390[1:(length(selic_4390)-1)], 10))
## [,1] [,2]
## [1,] 1.46 1.46
## [2,] 1.45 1.45
## [3,] 1.45 1.45
## [4,] 1.30 1.30
## [5,] 1.49 1.49
## [6,] 1.39 1.39
## [7,] 1.31 1.31
## [8,] 1.41 1.41
## [9,] 1.22 1.22
## [10,] 1.29 1.29
E acho que isso resolve uma parte do mistério. A série 4390 e a série Selic Overnight do site do IPEA de fato parecem ser a mesma série. E elas são o resultado da fórmula de juros compostos na série 11. Aisha: Agora resta saber por que a série 11 tem formato tão de escadinha.
Como é mais fácil de pegar os dados usando o pacote BETS, eu vou usar a série selic_4390 como a série “oficial”. Agora resta baixar os dados da série inteira (jul-94 até dez-17) e anualizar.
selic_4390 <- BETS.get("4390", from = "1994-07-01", to = "2017-12-31")
## Transformando em série anual
selic <- ((1+selic_4390/100)^(12)-1)*100
head(selic)
## Jul Aug Sep Oct Nov Dec
## 1994 121.95744 63.27210 56.99082 53.22268 61.40116 56.44736
tail(selic)
## Jul Aug Sep Oct Nov Dec
## 2017 10.033869 10.033869 7.956187 7.956187 7.058561 6.675963
length(selic)
## [1] 282
df2 <- data.frame(time(selic), as.numeric(selic))
names(df2) <- c("Data", "Selic")
p5 <- ggplot(df2, aes(Data, Selic)) +
geom_line() +
xlab("Data") +
ylab("Selic (%a.a.)") +
ggtitle("Taxa de Juros - Selic")
p5
Os valores do início da série estão meio ruins, mas depois eu me preocupo com isso se for o caso. Tem outra coisa me incomodando que os valores estão meio diferentes do que tem no site do Bacen: https://www.bcb.gov.br/Pec/Copom/Port/taxaSelic.asp. Aí não sei se eu não fiz merda na hora de converter a taxa.
Fim da digressão sobre a SELIC
As variáveis do modelo do Mumtaz são:
24364 que já tem ajuste sazonal (problema é que a série começa em janeiro de 2003 e acaba em outubro de 2017). Depois acabei pegando o PIB Mensal do Ipea Data deflacionado usando o IPCA de novembro de 2017. Essa deflação é feita com alguma bruxaria que não consegui identificar e por isso estou carregando o arquivo csv - não deu para baixar as séries do Bacen e deflacionar eu mesma e achar a mesma coisa.433 que já está em variação percentual mensal;7818 Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) (ela tem valores de 30/09/1999 a desembro de 2017). Depois comparei com a selic e como são praticamente iguais, fiquei com a Selic, mas fiz um teste com ela.3696 - taxa de câmbio livre - dólar americano (venda) fim de período mensal.Todas as variáveis exceto as de desigualdade, inflação e a taxa de juros estão em log differences (diferença dos logaritmos).
Bom, não tem PIB per capita, então tem que ser IBC-Br.
ibcbr <- BETS.get("24364", from = "2003-09-01", to = "2017-12-31")
df3 <- data.frame(time(ibcbr), as.numeric(ibcbr))
names(df3) <- c("Data", "IBCBr")
p6 <- ggplot(df3, aes(Data, IBCBr)) +
geom_line() +
xlab("Data") +
ylab("IBC-Br") +
ggtitle("Índice de Atividade Econômica do Bacen - com ajuste sazonal")
p6
LEMBRETE: Esse trem tem uma tendência monstra, depois tem que ver isso. Tirei a diferença da série em log.
LEMBRETE: A série começa em janeiro de 2003… Daí fica com 178 dados… tem que ver com o Guilherme depois.
Fuçando no site do Banco Central, eu achei uma série chamada PIB mensal - Valores correntes (R$ milhões) (4380).
Vou dar uma olhada nela:
PIBmensal <- BETS.get("4380", from = "1994-07-01", to = "2017-12-31")
data1 <- seq(as.Date("1994-07-01"), length = length(PIBmensal), by = "1 month")
df3 <- data.frame(data1, as.numeric(PIBmensal))
p6 <- ggplot(df3, aes(data1, PIBmensal)) +
geom_line() +
xlab("Data") +
ylab("PIB Mensal (R$ milhões)") +
ggtitle("PIB mensal - Valores correntes")
p6
## Don't know how to automatically pick scale for object of type ts. Defaulting to continuous.
# Como eu estou com preguiça de calcular o PIB deflacionado, peguei a série do site do IPEA pronta deflacionada pelo IPCA bonitinha
PIBmensal2 <- read.csv("C:\\Users\\Aishameriane\\OneDrive\\Documentos\\Mestrado Economia\\Teoria Macroeconômica II - 2017-2\\Artigo\\Aplicação\\PIBmensal.csv", header = T, sep = ";", dec = ",")
data1 <- seq(as.Date("1994-07-01"), length = length(PIBmensal2), by = "1 month")
df3 <- data.frame(data1, as.numeric(PIBmensal))
names(df3) <- c("Data", "PIBmensal")
p6 <- ggplot(df3, aes(Data, PIBmensal)) +
geom_line() +
xlab("Data") +
ylab("PIB Mensal (R$ milhões)") +
ggtitle("PIB mensal - Valores reais (base IPCA 11/2017)")
p6 + scale_y_continuous(labels = comma)
# Agora preciso transformar em índice
PIBmensal_in <- (PIBmensal/PIBmensal[1])*100
df3 <- data.frame(data1, as.numeric(PIBmensal_in))
names(df3) <- c("Data", "PIBmensal")
p6 <- ggplot(df3, aes(Data, PIBmensal)) +
geom_line() +
xlab("Data") +
ylab("PIB Mensal") +
ggtitle("PIB mensal índice - base = 2004/07")
p6
ibcbr <- BETS.get("24364", from = "2003-01-01", to = "2017-10-31")
data2 <- seq(as.Date("2003-01-01"), length = length(ibcbr), by = "1 month")
df4 <- data.frame(data2, as.numeric(ibcbr))
names(df4) <- c("Data", "IBCBr")
# Tentando juntar os pontos
test1 <- melt(data = df3, id.vars = "Data")
test2 <- melt(data = df4, id.vars = "Data")
test <- rbind(test1, test2)
ggplot(test, aes(Data, value, colour = variable)) +
geom_line(alpha = 1)+
labs(title="IBC-Br x PIB mensal", y = "", x= "Tempo", color = "Variável") +
scale_colour_brewer(palette = "Dark2") +
theme_bw()
PIBmensal_in <- ts(PIBmensal_in, start = c(1994, 07, 01), frequency = 12)
Update: Eu fiz alguns testes e o VAR com o PIB mensal “deu bom”. Como é importante saber fazer o que um aluno graduado em economia faz (MOURA, 2018), vou refazer as coisas em cima deflacionando eu mesma a série. Peguei um tutorial aqui. Update: Essa merda não funciona, vou continuar na ignorância mesmo, aff. Update: Havia um erro na merda da série do IPCA do IPEA data e aí os resultados entre o que eu calculava no Excel e o que tinha no Bacen davam diferentes. Eu e o Cássio mandamos um e-mail pro IPEA avisando. E agora eu consigo fazer tudo sem usar um csv. :D
Pra Aisha do futuro, eu tentei usar a fórmula e os valores deram diferentes do site do IPEA. Aí lendo as notas metodológicas, fala isso aqui: " Os índices IPCA, INPC e IGP-DI refletem o nível de preços médios de cada período, sendo inadequados para deflacionar séries referentes ao final do período. Visando minimizar tal problema, o sistema calcula esses índices de preço no final de cada período através do cálculo da média geométrica entre o índice do período e o índice do período seguinte." Calcular a média geométrica é tirar a raiz quadrada do produto, mas quando o índice dá negativo, dá ruim. Aí tentei fazer a média geométrica do IPCA acumulado, também não deu e eu desisti porque não tá sobrando tanto tempo assim. Agora deu.
PIBmensal <- BETS.get("4380", from = "1999-06-01", to = "2017-11-30")
data1 <- seq(as.Date("1994-07-01"), length = length(PIBmensal), by = "1 month")
df3 <- data.frame(data1, as.numeric(PIBmensal))
names(df3) <- c("Data", "PIBmensal")
p6 <- ggplot(df3, aes(Data, PIBmensal)) +
geom_line() +
xlab("Data") +
ylab("PIB Mensal (R$ milhões)") +
ggtitle("PIB mensal - Valores correntes")
p6
# O IPCA tem que ser pego desde 1993 caso queira confrontar com os valores do IPEA Data
ipca <- BETS.get("433", from = "1993-12-31", to = "2017-12-31")
ipca2 <- ipca
ipca2[1] <- 100
# Transformar em índice
for (i in 2:length(ipca2)){
ipca2[i] <- round(ipca2[(i-1)]*(1+ipca2[(i)]/100),2)
}
# Média geométrica
for (i in 1:(length(ipca2)-1)){
ipca2[i] <- sqrt(ipca2[i]*ipca2[(i+1)])
}
ipca2 <- ipca2[1:(length(ipca2)-1)]
# Mudança para base de nov/2017
for (i in 1:length(ipca2)) {
ipca2[i] <- tail(ipca2, n=1)/ipca2[i]
}
# Tirando as observações que não quero
ipca2 <- ts(ipca2, start = c(1993, 12, 31), frequency = 12)
ipca2 <- window(ipca2, c(1996, 6))
PIBmensal <- PIBmensal*ipca2
Para inflação usei a série 433, do IPCA. Tem duas coisas que podem ser feitas, uma é calcular o índice acumulado no ano (o que não parece interessante porque implica que todo início de ano teria uma queda) e a segunda é calcular o índice acumulado nos últimos 12 meses. Para fazer a segunda, tem que primeiro pegar todos os valores da série e dividir por 100 e depois acrescentar 1. Dessa nova série, são multiplicados os 12 valores anteriores (incluindo o atual), desconta um e multiplica por 100.
# Primeiro um gráfico simples
ipca <- BETS.get("433", from = "1994-07-01", to = "2017-12-31")
df4 <- data.frame(time(ipca), as.numeric(ipca))
names(df4) <- c("Data", "IPCA")
p7 <- ggplot(df4, aes(Data, IPCA)) +
geom_line() +
xlab("Data") +
ylab("IPCA") +
ggtitle("Inflação (IPCA) - mensal")
p7
# Agora a inflação acumulada em 12 meses
ipca_acum <- ipca/100 + 1
ipca_acum2 <- vector()
for (i in 12:282){
ipca_acum2[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}
ipca_acum2 <- ts(ipca_acum2, start = c(1995, 6, 1), frequency = 12)
df4 <- data.frame(time(ipca_acum2), as.numeric(ipca_acum2))
names(df4) <- c("Data", "IPCA")
p8 <- ggplot(df4, aes(Data, IPCA)) +
geom_line() +
xlab("Data") +
ylab("IPCA") +
ggtitle("Inflação (IPCA) - acumulada 12 meses")
p8
Eu conferi os valores com isso aqui e bateu. Já posso pedir a carteirinha de economista, yay! Posso nada, não consegui deflacionar o PIB. CONSEGUI, mwahahah.
Aisha: De novo tem aquelas valores ruins no início da série, acho que vou ter mesmo que eliminar uns valores.
Depois de falar com o Portela, cheguei nessa série aqui: 7818 Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) . Eu ainda não sei como que isso se relaciona a uma taxa de juros (no sentido de instrumento de política monetária), mas vou deixar os dados aqui guardados por enquanto.
swap90 <- BETS.get("7818", from = "1999-09-01", to = "2017-12-31")
df7 <- data.frame(time(swap90), as.numeric(swap90))
names(df7) <- c("Data", "DI90")
p10 <- ggplot(df7, aes(Data, DI90)) +
geom_line() +
xlab("Data") +
ylab("DI90 (%a.a.)") +
ggtitle("Taxa referencial de swaps DI pré 3 meses - fim de período")
p10
Vou tentar comparar com a SELIC. Será que essas duas taxas são comparáveis? Acho que sim porque eu passei SELIC para %a.a. também….
selic_4390 <- BETS.get("4390", from = "1999-09-01", to = "2017-12-31")
## Transformando em série anual
selic <- ((1+selic_4390/100)^(12)-1)*100
df10 <- data.frame(time(selic), as.numeric(selic))
names(df10) <- c("Data","Selic")
p11 <- ggplot(df10, aes(Data, Selic)) +
geom_line() +
xlab("Data") +
ylab("Selic (%a.a.)") +
ggtitle("Selic")
grid.arrange(p10, p11, ncol=1, nrow = 2)
df8 <- data.frame(df7, selic)
names(df8) <- c("Data", "DI90", "Selic")
df9 <- melt(data = df8, id.vars = "Data")
ggplot(df9, aes(Data, value, colour = variable)) +
geom_line(alpha = 1)+
labs(title="Taxas de juros", y = "", x= "Time", color = "Taxa") +
scale_colour_brewer(palette = "Dark2") +
theme_bw()
Elas são praticamente iguais, então vou ficar com a Selic. Update (18/01/2018) Guilherme mandou ficar com o swap DI. Update (22/01/2018) Agora que vi que essa coisa começa só em 99.
Do câmbio, tem essas séries no Bacen:
Como o Portela usa a segunda, vou usar a segunda (argumento de autoridade).
cambio <- BETS.get("3696", from = "1994-07-01", to = "2017-12-31")
df5 <- data.frame(time(cambio), as.numeric(cambio))
names(df5) <- c("Data", "Cambio")
p8 <- ggplot(df5, aes(Data, Cambio)) +
geom_line() +
xlab("Data") +
ylab("Tx. de Câmbio (u.m.c./US$)") +
ggtitle("Taxa de câmbio - Dólar americano (venda) - Fim de período - mensal")
p8
Realmente a partir de 99 que as coisas começaram a ficar diferentes.
Vou tentar calcular a taxa de crescimento usando a diferença dos logaritmos.
cambio_cres <- diff(log(cambio), 1)
df6 <- data.frame(time(cambio_cres), as.numeric(cambio_cres))
names(df6) <- c("Data", "Cambio")
p9 <- ggplot(df6, aes(Data, Cambio)) +
geom_line() +
xlab("Data") +
ylab("Variação") +
ggtitle("Taxa de variação da Taxa de câmbio - Dólar americano (venda) - Fim de período - mensal")
p9
Aí como estava com tempo sobrando, pensei em pegar o Ibovespa. Teoricamente teria alguma relação entre preço dos ativos e taxa de juros (esse é um dos canais de desigualdade).
symbols <- c("^BVSP")
getSymbols(symbols, src = 'yahoo', from = "1996-01-01",
auto.assign = TRUE, warnings = FALSE)
## 'getSymbols' currently uses auto.assign=TRUE by default, but will
## use auto.assign=FALSE in 0.5-0. You will still be able to use
## 'loadSymbols' to automatically load data. getOption("getSymbols.env")
## and getOption("getSymbols.auto.assign") will still be checked for
## alternate defaults.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.warning4.0"=FALSE). See ?getSymbols for details.
##
## WARNING: There have been significant changes to Yahoo Finance data.
## Please see the Warning section of '?getSymbols.yahoo' for details.
##
## This message is shown once per session and may be disabled by setting
## options("getSymbols.yahoo.warning"=FALSE).
## Warning: ^BVSP contains missing values. Some functions will not work if
## objects contain missing values in the middle of the series. Consider using
## na.omit(), na.approx(), na.fill(), etc to remove or replace them.
## [1] "BVSP"
prices <- merge(BVSP)
prices <- prices[,grepl( "Close" , names(prices) )]
colnames(prices) <- symbols
returns <- na.omit(monthlyReturn(prices, type = "log"))
## Warning in to_period(xx, period = on.opts[[period]], ...): missing values
## removed from data
plot.xts(returns)
ibovespa <- ts(returns, start = c(1996,1,1), frequency = 12)
df1 <- data.frame(seq(as.Date("1996-01-01"), length = length(ibovespa), by = "1 month"), ibovespa)
names(df1) <- c("Data", "Ibovespa")
p9 <- ggplot(df1, aes(Data, Ibovespa)) +
geom_line() +
xlab("Data") +
ylab("Ibovespa") +
scale_x_date(date_breaks = "12 months")+
ggtitle("Log retornos mensais do Ibovespa")
p9 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
A lista de dados que tenho até agora é:
Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do trabalho (7620)Receitas tributárias - Regime de competência - Imposto de renda - Retido na fonte - Rendimento do capital (7621) Não foi feita nenhuma transformação nos dados.Taxa de juros - Selic acumulada no mês (4390) Eu comparei com a taxa referencial de swaps DI pré-fixada e ambas são praticamente a mesma coisa, então fiquei com a Selic mesmo. A única transformação feita (até agora) foi passar de taxa mensal para anual usando a fórmula: \(\left((1+tx/100)^12 -1\right)*100\).
Taxa referencial de swaps DI pré-fixada (BM&F) - Prazo de 90 dias (fim de período) (7818) Sem transformação (a taxa já é ao ano). O problema é que essa série começa apenas em 1999.
Índice nacional de preços ao consumidor-amplo (IPCA) (433) Como ele está em variação percentual mensal, foi necessário transformar para acumulado dos últimos 12 meses utilizando a fórmula: \[IPCA_i = \left[\left(\prod\limits_{j=i-11}^i \left(\frac{IPCA_j}{100}+1\right) \right) -1\right]*100\]Índice de Atividade Econômica do Banco Central (IBC-Br) - com ajuste sazonal (24364) Taxa de câmbio - Livre - Dólar americano (venda) - Fim de período - mensal (3696) Como o Mumtaz disse que usou log diferença, eu calculei o log e tirei primeira diferença nessa série. rm(list = ls())
# Auxiliary variables, so I don't need to bother when something changes
inicio <- "2003-02-01"
fim <- "2017-10-31"
# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries
trabalho <- BETS.get("7620", from = inicio, to = fim)
capital <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho
# Usando ndiffs(ibcbr,test="adf",alpha = 0.1) se encontra que o ibcbr é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
ibcbr_raw <- BETS.get("24364", from = inicio_cambio, to = fim)
ibcbr <- diff(log(ibcbr_raw), 1)
cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio <- diff(log(cambio_raw), 1)
selic_4390 <- BETS.get("4390", from = inicio, to = fim)
selic <- ((1+selic_4390/100)^(12)-1)*100
ipca_raw <- BETS.get("433", from = inicio_ipca, to = fim)
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()
final <- length(ipca_acum)
for (i in 12:final){
ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}
ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))
ipca <- ts(ipca, start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD
# Plotting some nice graphs
df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho, selic, ibcbr, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "IBCBr", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")
ggplot(df2, aes(Data, value, colour = variable)) +
geom_line(alpha = 1)+
labs(title="Variáveis do modelo", y = "", x= "Time", color = "Variável") +
scale_colour_brewer(palette = "Dark2") +
theme_bw()
cores <- brewer.pal(5, "Dark2")
# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Capital trabalho") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p2 <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[2])+
scale_y_continuous(name="Selic") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p3 <- ggplot(df2[which(df2$variable == "IBCBr"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[3])+
scale_y_continuous(name="IBC-Br") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[4])+
scale_y_continuous(name="Tx. Câmbio (var)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p5 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[5])+
scale_y_continuous(name="IPCA") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
theme_bw()
p5 <- p5 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5)
Acrescentando o PIB e o Swap DI 3 meses:
rm(list = ls())
# Auxiliary variables, so I don't need to bother when something changes
inicio <- "1996-01-01"
fim <- "2017-11-30"
inicio_deflator <- "1993-12-31"
# Don't mess with this code
inicio_cambio <- paste(seq(as.Date(inicio), length = 2, by = "-1 month")[2]) # 1 mês antes do início das outras séries
inicio_ipca2 <- paste(seq(as.Date(inicio), length = 12, by = "-1 month")[12]) # 12 meses antes do início das outras séries
fim_ipca <- paste(seq(as.Date(fim), length = 2, by = "1 month")[2]) # 1 mês depois do fim das outras séries
trabalho <- BETS.get("7620", from = inicio, to = fim)
capital <- BETS.get("7621", from = inicio, to = fim)
capital_trabalho <- capital/trabalho
ipca_raw <- BETS.get("433", from = inicio_deflator, to = fim_ipca)
ipca2 <- ipca_raw
ipca2[1] <- 100
# Transformar em índice
for (i in 2:length(ipca2)){
ipca2[i] <- round(ipca2[(i-1)]*(1+ipca2[(i)]/100),2)
}
# Média geométrica
for (i in 1:(length(ipca2)-1)){
ipca2[i] <- sqrt(ipca2[i]*ipca2[(i+1)])
}
# Mudança para base de nov/2017
for (i in 1:length(ipca2)) {
ipca2[i] <- tail(ipca2, n=2)[1]/ipca2[i]
}
# Tirando as observações que não quero
ipca2 <- ts(ipca2, start = c(1993, 12, 31), frequency = 12)
ano_ipca <- as.numeric(substr(inicio_cambio, start = 1, stop = 4))
mes_ipca <- as.numeric(substr(inicio_cambio, start = 6, stop = 7))
ipca2 <- window(ipca2, c(ano_ipca, mes_ipca))
ano_ipca <- as.numeric(substr(inicio_ipca2, start = 1, stop = 4))
mes_ipca <- as.numeric(substr(inicio_ipca2, start = 6, stop = 7))
ipca_raw <- window(ipca_raw, c(ano_ipca, mes_ipca))
# Tira dezembro 2017
ipca_raw <- ipca_raw[1:(length(ipca_raw)-1)]
ipca_acum <- ipca_raw/100 + 1
ipca <- vector()
final <- length(ipca_acum)
for (i in 12:final){
ipca[(i-11)] <- (prod(ipca_acum[(i-11):i])-1)*100
}
ano <- as.numeric(substr(inicio, start = 1, stop = 4))
mes <- as.numeric(substr(inicio, start = 6, stop = 7))
dia <- as.numeric(substr(inicio, start = 9, stop = 10))
ipca <- ts(ipca, start = c(ano, mes, dia), frequency = 12) # Date format YYYY MM DD
# Usando ndiffs(PIBmensal,test="adf",alpha = 0.05) se encontra que o PIBmensal é não estacionário.
# Como todas as outras séries são estacionárias, eu vou fazer também a diferença do log, igual no câmbio
PIBmensal <- BETS.get("4380", from = inicio_cambio, to = fim)
PIBmensal <- PIBmensal * ipca2
PIB <- diff(log(PIBmensal), 1)
cambio_raw <- BETS.get("3696", from = inicio_cambio, to = fim)
cambio <- diff(log(cambio_raw), 1)
swap90 <- BETS.get("7818", from = inicio, to = fim)
swap90 <- ts.union(swap90, ipca)
swap90 <- swap90[,1]
selic_4390 <- BETS.get("4390", from = inicio, to = fim)
selic <- ((1+selic_4390/100)^(12)-1)*100
# Plotting some nice graphs
df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho, swap90, PIB, cambio, selic, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Swap90", "PIB", "Cambio", "selic", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")
ggplot(df2, aes(Data, value, colour = variable)) +
geom_line(alpha = 1)+
labs(title="Variáveis do modelo", y = "", x= "Time", color = "Variável") +
scale_colour_brewer(palette = "Dark2") +
theme_bw()
cores <- brewer.pal(6, "Dark2")
# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Capital trabalho") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p2 <- ggplot(df2[which(df2$variable == "Swap90"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[2])+
scale_y_continuous(name="SwapDI90 (%a.a.)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p3 <- ggplot(df2[which(df2$variable == "PIB"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[3])+
scale_y_continuous(name="PIB mensal (tx.cresc)") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p3 <- p3 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[4])+
scale_y_continuous(name="Tx. Câmbio (var)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p5 <- ggplot(df2[which(df2$variable == "selic"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[5])+
scale_y_continuous(name="Selic (%a.a)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p5 <- p5 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
p6 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[6])+
scale_y_continuous(name="IPCA") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
theme_bw()
p6 <- p6 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 8))
grid.arrange(p1, p2, p4, p5, p3, p6, ncol=2, nrow = 3)
grid.arrange(p1, p4, p3, ncol=1, nrow = 3)
grid.arrange(p2, p5, p6, ncol = 1, nrow = 3)
Agora começa outro problema para tratar os dados: para a razão capital trabalho original, tem 30 observações de antes do plano real que tem valores muito altos. Depois de eliminar elas, tem duas séries com tendência e sazonalidade. A razão entre elas não aparenta ter tendência (fiz um teste que indicou que a série precisaria de uma diferença para ficar estacionária), porém parece ter uma componente sazonal persistente.
Resolvi tentar duas coisas: a primeira foi fazer manualmente a remoção de tendência e sazonalidade e a segunda foi usar um filtro pronto. Meu filtro manual ficou uma merda, então vou ficar com o automático mesmo.
Vou deixar aqui e qualquer coisa apago no final. Os chunks não estão sendo rodados. Usei este tutorial - https://anomaly.io/seasonal-trend-decomposition-in-r/.
Não sei se precisa. Como os dados são mensais, vou usar uma janela móvel de tamanho 12.
# Checa número de diferenciações necessárias para tornar séries estacionárias
ndiffs(capital_trabalho,test="adf",alpha = 0.1)
# De acordo com o teste, seria necessário diferenciar uma vez a série para que ela fique estacionária
# Remove a tendência ajustando uma janela móvel
trend_capital_trabalho = ma(capital_trabalho, order = 12, centre = T)
plot(as.ts(capital_trabalho))
lines(trend_capital_trabalho)
plot(as.ts(trend_capital_trabalho))
df <- as.data.frame(cbind(time(trend_capital_trabalho),trend_capital_trabalho, capital_trabalho))
head(df)
names(df) <- c("Data", "Tendencia", "CapitalTrabalho")
ggplot(df, aes(Data)) +
geom_line(aes(y = CapitalTrabalho, colour = "Capital/Trabalho")) +
geom_line(aes(y = Tendencia, colour = "Tendência"))
m_capital_trabalho = t(matrix(data = capital_trabalho, nrow = 12))
seasonal_capital_trabalho = colMeans(m_capital_trabalho, na.rm = T)
plot(as.ts(rep(seasonal_capital_trabalho,12)))
Vamos calcular o “random noise” usando \(Resíduo = Série - Sazonalidade - Tendência\).
residuo = capital_trabalho - seasonal_capital_trabalho - trend_capital_trabalho
plot(as.ts(residuo))
Ainda assim está meio estranho… Vou tentar usar um filtro automático.
O X-13ARIMA-SEATS é um algoritmo de ajuste sazonal desenvolvido pelo United States Census Bureau. Sua implementação no R é feita pelo pacote seasonal. Um manual pode ser visto aqui: https://cran.r-project.org/web/packages/seasonal/vignettes/seas.pdf. Usei muito este documento aqui: file:///C:/Users/Aishameriane/Desktop/Ferreira_Mattos_2016.pdf e também o teste de quebra estrutural que tem aqui https://pythonandr.com/2016/11/08/endogenously-detecting-structural-breaks-in-a-time-series-implementation-in-r/ .
Update (18/01/2018) O Guilherme sugeriu passar um filtro na quebra da razão capital trabalho em 1999. Primeiro vou tirar a sazonalidade da série toda, ver como fica o gráfico, quebrar em dois períodos a série original, depois ajudar o filtro sazonal para cada uma das séries.
#######################
# Capital trabalho #
#######################
cores <- brewer.pal(6, "Dark2")
# Fazendo um teste incluindo os dados desde 1996 (que depois quero fazer o VAR com isso)
trabalho <- BETS.get("7620", from = "1996-01-01", to = "2017-11-01")
capital <- BETS.get("7621", from = "1996-01-01", to = "2017-11-01")
capital_trabalho <- capital/trabalho
m <- seas(x = capital_trabalho, transform.function = "none", regression.aictest = NULL)
qs(m)
## qs p-val
## qsori 285.61105 0.00000
## qsorievadj 378.83903 0.00000
## qsrsd 0.00000 1.00000
## qssadj 46.77278 0.00000
## qssadjevadj 0.00000 1.00000
## qsirr 48.12373 0.00000
## qsirrevadj 0.00000 1.00000
## qssori 127.09650 0.00000
## qssorievadj 146.78226 0.00000
## qssrsd 0.00000 1.00000
## qsssadj 11.73690 0.00283
## qsssadjevadj 0.00000 1.00000
## qssirr 9.09815 0.01058
## qssirrevadj 0.00000 1.00000
plot(m)
summary(m)
##
## Call:
## seas(x = capital_trabalho, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Constant 0.023146 0.005041 4.592 4.40e-06 ***
## AO1996.Jan -0.657991 0.085330 -7.711 1.25e-14 ***
## LS1996.Jul -0.215914 0.035320 -6.113 9.77e-10 ***
## LS1998.Jan 0.851747 0.060992 13.965 < 2e-16 ***
## LS1998.Feb -0.456248 0.075658 -6.030 1.64e-09 ***
## LS1998.Apr -0.343722 0.057834 -5.943 2.79e-09 ***
## LS1998.Jul 0.339913 0.036980 9.192 < 2e-16 ***
## AO1998.Aug 0.373609 0.060735 6.151 7.68e-10 ***
## AO1999.Feb 0.712392 0.065041 10.953 < 2e-16 ***
## AO1999.Mar 0.518619 0.065041 7.974 1.54e-15 ***
## LS1999.Dec -0.237452 0.024535 -9.678 < 2e-16 ***
## AO2001.Jun 0.297887 0.059118 5.039 4.68e-07 ***
## AO2001.Sep 0.267972 0.059117 4.533 5.82e-06 ***
## AO2001.Oct 0.301303 0.065105 4.628 3.69e-06 ***
## AO2002.Oct 0.580093 0.065104 8.910 < 2e-16 ***
## AO2003.Feb 0.330971 0.059131 5.597 2.18e-08 ***
## AO2003.Dec -0.271940 0.060961 -4.461 8.16e-06 ***
## LS2004.Mar -0.196492 0.025369 -7.745 9.53e-15 ***
## LS2004.Nov -0.182872 0.026266 -6.962 3.35e-12 ***
## AO2005.Jun 1.001053 0.072560 13.796 < 2e-16 ***
## AO2005.Dec 0.394815 0.059238 6.665 2.65e-11 ***
## AO2006.Jan 0.320392 0.065114 4.920 8.63e-07 ***
## AO2006.Jun 0.996463 0.081624 12.208 < 2e-16 ***
## AO2007.Jan 0.265836 0.065064 4.086 4.39e-05 ***
## AO2007.Jun 0.691402 0.080924 8.544 < 2e-16 ***
## AO2008.Jun 0.244229 0.070170 3.481 0.00050 ***
## LS2009.Aug -0.099893 0.023161 -4.313 1.61e-05 ***
## AO2011.Oct 0.261397 0.059084 4.424 9.68e-06 ***
## AO2011.Dec 0.241221 0.059084 4.083 4.45e-05 ***
## AO2013.Jun -0.411605 0.065042 -6.328 2.48e-10 ***
## AO2014.Jun -0.335137 0.065042 -5.153 2.57e-07 ***
## MA-Seasonal-12 0.163729 0.063369 2.584 0.00977 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (0 0 0)(0 1 1) Obs.: 263 Transform: none
## AICc: -495.2, BIC: -389.2 QS (no seasonality in final):46.77 ***
## Box-Ljung (no autocorr.): 38.53 * Shapiro (normality): 0.9772 ***
## Messages generated by X-13:
## Notes:
## - Unable to test AO1998.Jan due to regression matrix
## singularity.
sliding <- series(m, "sfs")
## specs have been added to the model: slidingspans
sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100
## [1] 0
ggmonthplot(capital_trabalho) +
theme_bw() +
scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "Nenhuma transformação")
# Gráfico por ano e mês
ggseasonplot(capital_trabalho, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
labs(title = "Média mensal da razão capital-trabalho", x = "Mês", subtitle = "Nenhuma transformação")
# Usa a interface gráfica
#view(m)
# Salva a série final em uma nova variável
capital_trabalho2 <- final(m)
#plot(capital_trabalho2)
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho2), by = "1 month"), capital_trabalho, capital_trabalho2)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = capital_trabalho, colour = "Original")) +
geom_line(aes(y = capital_trabalho2, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
scale_y_continuous(name = "Capital/Trabalho")+
labs(color="Variável") +
theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1
# Dividindo a série em pedaços distintos:
# pega os pontos de quebras
bp_ts <- breakpoints(capital_trabalho ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = capital_trabalho ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 97
## m = 2 39 97
## m = 3 39 97 221
## m = 4 39 97 138 221
## m = 5 39 97 138 185 224
##
## Corresponding to breakdates:
##
## m = 1 2004(1)
## m = 2 1999(3) 2004(1)
## m = 3 1999(3) 2004(1) 2014(5)
## m = 4 1999(3) 2004(1) 2007(6) 2014(5)
## m = 5 1999(3) 2004(1) 2007(6) 2011(5) 2014(8)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 22.23 20.07 19.79 19.63 19.42 19.40
## BIC 107.76 92.01 99.40 108.37 116.75 127.56
# Calcula o intervalo de confiança
ci_ts <- confint(bp_ts)
# Monta o dataframe
Data <- seq(as.Date("1996-01-01"), length = length(capital_trabalho), by = "1 month")
df1 <- melt(data.frame(Data,capital_trabalho), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")
p1 <- ggplot(df1[which(df1$Variavel == "capital_trabalho"),], aes(Data, Valor, colour = Variavel)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Razão Capital/Trabalho") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(capital_trabalho), xend = Data[ci_ts$confint[3]], yend = min(capital_trabalho)))+
theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1
# Salvando duas séries para passar o seasonal
capital_trabalho_a <- capital_trabalho[1:bp_ts$breakpoints]
capital_trabalho_b <- capital_trabalho[(bp_ts$breakpoints+1):length(capital_trabalho)]
capital_trabalho_a <- ts(capital_trabalho_a, start = c(1996, 01, 01), frequency = 12)
capital_trabalho_b <- ts(capital_trabalho_b, start = c(2004, 02, 01), frequency = 12)
m <- seas(x = capital_trabalho_a, transform.function = "none", regression.aictest = NULL)
qs(m)
## qs p-val
## qsori 34.12733 0.00000
## qsorievadj 91.53477 0.00000
## qsrsd 0.00000 1.00000
## qssadj 0.00000 1.00000
## qssadjevadj 0.00000 1.00000
## qsirr 0.00078 0.99961
## qsirrevadj 0.00000 1.00000
## qssori 33.13775 0.00000
## qssorievadj 85.86391 0.00000
## qssrsd 0.00000 1.00000
## qsssadj 0.00000 1.00000
## qsssadjevadj 0.00000 1.00000
## qssirr 0.05069 0.97497
## qssirrevadj 0.00000 1.00000
plot(m)
summary(m)
##
## Call:
## seas(x = capital_trabalho_a, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AO1996.Jan -0.35561 0.10377 -3.427 0.000611 ***
## AO1998.Jan 0.67046 0.08634 7.766 8.12e-15 ***
## LS1998.Jul 0.28775 0.05339 5.390 7.06e-08 ***
## AO1998.Aug 0.29547 0.08564 3.450 0.000561 ***
## AO1999.Feb 0.68086 0.09704 7.016 2.28e-12 ***
## AO1999.Mar 0.37531 0.09644 3.892 9.96e-05 ***
## AO2002.Oct 0.49217 0.08564 5.747 9.10e-09 ***
## AR-Nonseasonal-01 0.64477 0.07100 9.081 < 2e-16 ***
## MA-Seasonal-12 0.99965 0.08975 11.138 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (1 0 0)(0 1 1) Obs.: 97 Transform: none
## AICc: -109.8, BIC: -88.36 QS (no seasonality in final): 0
## Box-Ljung (no autocorr.): 20.02 Shapiro (normality): 0.9765 .
sliding <- series(m, "sfs")
## specs have been added to the model: slidingspans
sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100
## Warning in is.na(object): is.na() aplicado a um objeto diferente de lista
## ou vetor de tipo 'NULL'
## [1] NaN
#view(m)
capital_trabalho_2a <- final(m)
ggmonthplot(capital_trabalho_2a) +
theme_bw() +
scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "De Jan/96 a Jan/04")
# Gráfico por ano e mês
ggseasonplot(capital_trabalho_2a, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
labs(title = "Média mensal da razão capital-trabalho", x = "Mês", subtitle = "De Jan/96 a Jan/04")
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho_a), by = "1 month"), capital_trabalho_a, capital_trabalho_2a)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = capital_trabalho_a, colour = "Original")) +
geom_line(aes(y = capital_trabalho_2a, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Capital/Trabalho")+
labs(color="Variável") +
theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1
# Dividindo a série em pedaços distintos:
# pega os pontos de quebras
bp_ts <- breakpoints(capital_trabalho_a ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = capital_trabalho_a ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 24
## m = 2 24 42
## m = 3 24 49 65
## m = 4 24 49 65 81
## m = 5 24 39 53 67 81
##
## Corresponding to breakdates:
##
## m = 1 1997(12)
## m = 2 1997(12) 1999(6)
## m = 3 1997(12) 2000(1) 2001(5)
## m = 4 1997(12) 2000(1) 2001(5) 2002(9)
## m = 5 1997(12) 1999(3) 2000(5) 2001(7) 2002(9)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 6.505 4.558 4.078 3.834 3.791 3.990
## BIC 22.318 -3.045 -4.687 -1.523 6.530 20.643
# Calcula o intervalo de confiança
ci_ts <- confint(bp_ts)
#############
m <- seas(x = capital_trabalho_b, transform.function = "none", regression.aictest = NULL)
qs(m)
## qs p-val
## qsori 229.41027 0.00000
## qsorievadj 277.61642 0.00000
## qsrsd 0.00000 1.00000
## qssadj 3.02991 0.21982
## qssadjevadj 0.00000 1.00000
## qsirr 2.37126 0.30555
## qsirrevadj 0.00000 1.00000
## qssori 124.60633 0.00000
## qssorievadj 149.00319 0.00000
## qssrsd 0.00000 1.00000
## qsssadj 6.03848 0.04884
## qsssadjevadj 0.00000 1.00000
## qssirr 7.36388 0.02517
## qssirrevadj 0.00000 1.00000
plot(m)
summary(m)
##
## Call:
## seas(x = capital_trabalho_b, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AO2004.Jun -0.75270 0.08042 -9.359 < 2e-16 ***
## LS2004.Nov -0.33098 0.04195 -7.889 3.03e-15 ***
## AO2004.Dec -0.39821 0.05723 -6.959 3.44e-12 ***
## LS2005.Mar 0.15205 0.03782 4.020 5.82e-05 ***
## AO2005.Jun 0.37303 0.07240 5.152 2.57e-07 ***
## AO2006.Jun 0.51048 0.06305 8.097 5.63e-16 ***
## AO2006.Dec -0.22518 0.04357 -5.168 2.36e-07 ***
## AO2007.Jun 0.35674 0.05244 6.803 1.02e-11 ***
## AO2008.Jan -0.19683 0.04278 -4.601 4.21e-06 ***
## AO2010.Jun -0.19763 0.04285 -4.612 3.98e-06 ***
## AO2011.Oct 0.23059 0.04291 5.374 7.69e-08 ***
## AO2011.Dec 0.22011 0.04291 5.130 2.90e-07 ***
## AO2013.Jun -0.40288 0.04538 -8.878 < 2e-16 ***
## AO2014.Jun -0.31472 0.04536 -6.939 3.96e-12 ***
## AO2016.Dec 0.20542 0.05233 3.925 8.67e-05 ***
## MA-Nonseasonal-01 0.73304 0.05551 13.204 < 2e-16 ***
## MA-Seasonal-12 0.33591 0.07718 4.352 1.35e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (0 1 1)(0 1 1) Obs.: 166 Transform: none
## AICc: -403.3, BIC: -353.9 QS (no seasonality in final): 3.03
## Box-Ljung (no autocorr.): 12.88 Shapiro (normality): 0.9618 ***
sliding <- series(m, "sfs")
## specs have been added to the model: slidingspans
sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100
## [1] 0
#view(m)
capital_trabalho_2b <- final(m)
ggmonthplot(capital_trabalho_2b) +
theme_bw() +
scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
labs(title = "Média diária e mensal da razão capital-trabalho", x = "Mês", subtitle = "De Fev/04 a Nov/17")
# Gráfico por ano e mês
ggseasonplot(capital_trabalho_2b, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Razão Capital/Trabalho (% - média)") +
labs(title = "Média mensal da razão capital-trabalho", x = "Mês", subtitle = "De Fev/04 a Nov/17")
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("2004-02-01"), length = length(capital_trabalho_b), by = "1 month"), capital_trabalho_b, capital_trabalho_2b)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = capital_trabalho_b, colour = "Original")) +
geom_line(aes(y = capital_trabalho_2b, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Capital/Trabalho")+
labs(color="Variável") +
theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1
# Dividindo a série em pedaços distintos:
# pega os pontos de quebras
bp_ts <- breakpoints(capital_trabalho_b ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = capital_trabalho_b ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 136
## m = 2 29 136
## m = 3 29 65 130
## m = 4 29 60 88 136
## m = 5 29 60 88 112 136
##
## Corresponding to breakdates:
##
## m = 1 2015(5)
## m = 2 2006(6) 2015(5)
## m = 3 2006(6) 2009(6) 2014(11)
## m = 4 2006(6) 2009(1) 2011(5) 2015(5)
## m = 5 2006(6) 2009(1) 2011(5) 2013(5) 2015(5)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 13.57 13.33 13.12 13.09 13.02 13.02
## BIC 65.60 72.89 80.47 90.37 99.64 109.83
# Calcula o intervalo de confiança
#ci_ts <- confint(bp_ts)
###############
# Juntando novamente as duas séries
###############
comb <- ts.union(capital_trabalho_2a, capital_trabalho_2b)
capital_trabalho_final <- pmin(comb[,1], comb[,2], na.rm = TRUE)
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho_final), by = "1 month"), capital_trabalho, capital_trabalho_final)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = capital_trabalho, colour = "Original")) +
geom_line(aes(y = capital_trabalho_final, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Razão Capital/Trabalho")+
labs(color="Variável") +
theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho2), by = "1 month"), capital_trabalho, capital_trabalho2, capital_trabalho_final)
names(df) <- c("Data", "capital_trabalho", "capital_trabalho2", "capital_trabalho_final")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = capital_trabalho, colour = "Original")) +
geom_line(aes(y = capital_trabalho2, colour = "Automático")) +
geom_line(aes(y = capital_trabalho_final, colour = "Manual")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Razão Capital/Trabalho")+
labs(color="Método") +
theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1
#######################
# Selic #
#######################
cores <- brewer.pal(6, "Dark2")
selic_4390 <- BETS.get("4390", from = "1996-01-01", to = "2017-11-30")
## Transformando em série anual
selic <- ((1+selic_4390/100)^(12)-1)*100
ggmonthplot(selic) +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a - média)") +
labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")
# Gráfico por ano e mês
ggseasonplot(selic, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Nenhuma transformação")
m <- seas(x = selic, transform.function = "none", regression.aictest = NULL)
#plot(m)
summary(m)
##
## Call:
## seas(x = selic, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## LS1996.Mar -4.09295 0.97485 -4.199 2.69e-05 ***
## LS1997.Nov 21.72884 0.95722 22.700 < 2e-16 ***
## LS1998.Jan -6.54038 0.99031 -6.604 3.99e-11 ***
## LS1998.Feb -5.49776 0.98955 -5.556 2.76e-08 ***
## LS1998.Apr -7.40308 0.95723 -7.734 1.04e-14 ***
## LS1998.Sep 15.32547 1.03446 14.815 < 2e-16 ***
## AO1998.Oct 7.50672 0.97560 7.694 1.42e-14 ***
## AO1998.Nov 3.92048 0.85822 4.568 4.92e-06 ***
## AO1999.Jan -4.09824 0.92123 -4.449 8.64e-06 ***
## AO1999.Mar 14.93677 0.85625 17.444 < 2e-16 ***
## LS1999.May -5.36681 1.05880 -5.069 4.00e-07 ***
## LS1999.Jun -4.48989 0.92662 -4.845 1.26e-06 ***
## AR-Nonseasonal-01 -0.25889 0.10018 -2.584 0.00976 **
## AR-Nonseasonal-02 0.13521 0.06290 2.150 0.03159 *
## AR-Nonseasonal-03 0.49771 0.05109 9.742 < 2e-16 ***
## MA-Nonseasonal-01 0.14349 0.11460 1.252 0.21053
## MA-Seasonal-12 1.00000 0.04180 23.923 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (3 1 1)(0 1 1) Obs.: 263 Transform: none
## AICc: 799.2, BIC: 859.7 QS (no seasonality in final):2.135
## Box-Ljung (no autocorr.): 59.74 *** Shapiro (normality): 0.9856 **
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
## been found in one or more of the estimated spectra.
##
## Notes:
## - Unable to test LS1998.Jan due to regression matrix
## singularity.
## - Unable to test LS1999.Jun due to regression matrix
## singularity.
## - Unable to test LS1998.Feb due to regression matrix
## singularity.
## - Unable to test LS1998.Feb due to regression matrix
## singularity.
## - Unable to test AO1998.Jan due to regression matrix
## singularity.
## - Unable to test LS1999.May due to regression matrix
## singularity.
qs(m)
## qs p-val
## qsori 23.76401 0.00001
## qsorievadj 96.16089 0.00000
## qsrsd 0.00000 1.00000
## qssadj 2.13450 0.34395
## qssadjevadj 18.30242 0.00011
## qsirr 0.00000 1.00000
## qsirrevadj 0.00000 1.00000
## qssori 33.30455 0.00000
## qssorievadj 33.30455 0.00000
## qssrsd 0.26400 0.87634
## qsssadj 6.61682 0.03657
## qsssadjevadj 6.61682 0.03657
## qssirr 0.00000 1.00000
## qssirrevadj 0.00000 1.00000
#summary(m) diz Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks se for a série a partir de 99
# Usa a interface gráfica
#view(m)
sliding <- series(m, "sfs")
## specs have been added to the model: slidingspans
sum(sliding[,"Max_._DIFF"] > 3, na.rm = T)/length(na.omit(sliding[,"Max_._DIFF"])) * 100
## [1] 0
# Salva a série final em uma nova variável
selic2 <- final(m)
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date(inicio), length = length(selic2), by = "1 month"), selic, selic2)
names(df) <- c("Data", "Original", "Série Filtrada")
p2 <- ggplot(df, aes(Data)) +
geom_line(aes(y = selic, colour = "Original")) +
geom_line(aes(y = selic2, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Selic (%a.a.)")+
labs(color="Variável") +
theme_bw()+theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p2
# Dividindo a série em pedaços distintos:
# pega os pontos de quebras
bp_ts <- breakpoints(selic ~ 1)
# o valor com o menor BIC é o escolhido
summary(bp_ts)
##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = selic ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 121
## m = 2 41 128
## m = 3 41 128 222
## m = 4 41 123 162 222
## m = 5 41 81 121 160 222
##
## Corresponding to breakdates:
##
## m = 1 2006(1)
## m = 2 1999(5) 2006(8)
## m = 3 1999(5) 2006(8) 2014(6)
## m = 4 1999(5) 2006(3) 2009(6) 2014(6)
## m = 5 1999(5) 2002(9) 2006(1) 2009(4) 2014(6)
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 13285 5841 3529 3417 3262 3211
## BIC 1789 1584 1463 1465 1464 1471
# Calcula o intervalo de confiança
ci_ts <- confint(bp_ts)
Data <- seq(as.Date("1996-01-01"), length = length(selic), by = "1 month")
df1 <- melt(data.frame(Data,selic), id.vars = "Data")
names(df1) <- c("Data", "Variavel", "Valor")
p1 <- ggplot(df1[which(df1$Variavel == "selic"),], aes(Data, Valor, colour = Variavel)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Selic (%a.a.)") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
geom_vline(xintercept = Data[bp_ts$breakpoints], linetype="longdash")+
geom_segment(aes(x = Data[ci_ts$confint[1]], y = min(selic), xend = Data[ci_ts$confint[5]], yend = min(selic)))+
geom_segment(aes(x = Data[ci_ts$confint[2]], y = min(selic), xend = Data[ci_ts$confint[6]], yend = min(selic)))+
theme_bw()
p1 <- p1 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), legend.position = "none")
p1
###########
# Dividindo a selic em três pedaços
##########
# Salvando três séries para passar o seasonal
selic_a <- selic[1:(bp_ts$breakpoints[1]-1)]
selic_b <- selic[(bp_ts$breakpoints[1]):(bp_ts$breakpoints[2]-1)]
selic_c <- selic[(bp_ts$breakpoints[2]):(length(selic))]
selic_a <- ts(selic_a, start = c(1996, 01, 01), frequency = 12)
selic_b <- ts(selic_b, start = c(1999, 05, 01), frequency = 12)
selic_c <- ts(selic_c, start = c(2006, 08, 01), frequency = 12)
#######
# Primeiro pedaço
#########
ggmonthplot(selic_a) +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a - média)") +
labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Jan/96 a Abr/99")
# Gráfico por ano e mês
ggseasonplot(selic_a, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Jan/96 a Abr/99")
m <- seas(x = selic_a, transform.function = "none", regression.aictest = NULL)
plot(m)
summary(m)
##
## Call:
## seas(x = selic_a, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Constant -0.6762 0.3314 -2.040 0.0413 *
## LS1997.Nov 21.4637 1.3511 15.886 < 2e-16 ***
## LS1998.Feb -9.0236 1.2976 -6.954 3.55e-12 ***
## LS1998.Apr -7.0057 1.3511 -5.185 2.16e-07 ***
## LS1998.Sep 15.4392 1.4204 10.870 < 2e-16 ***
## AO1998.Oct 5.8520 1.0364 5.646 1.64e-08 ***
## AO1999.Mar 15.4913 1.0364 14.947 < 2e-16 ***
## AR-Seasonal-12 0.4275 0.1998 2.140 0.0323 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (0 1 0)(1 0 0) Obs.: 40 Transform: none
## AICc: 164, BIC: 172.7 QS (no seasonality in final): 0
## Box-Ljung (no autocorr.): 31.97 Shapiro (normality): 0.9405 *
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
## because the spectrum of the prior adjusted series (Table
## B1) has no visually significant seasonal peaks.
##
## Notes:
## - Insufficient data to compute average forecast error
## diagnostic.
## - The program cannot perform hypothesis tests for kurtosis
## on less than 50 observations.
qs(m)
## qs p-val
## qsori 0.02853 0.98584
## qsorievadj 4.23955 0.12006
## qsrsd 1.32801 0.51479
## qssadj 0.00000 1.00000
## qssadjevadj 0.00000 1.00000
## qsirr 0.01122 0.99441
## qsirrevadj 0.00000 1.00000
#summary(m) #diz Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks
# Deixei ela como está
selic_2a <- selic_a
#######
# Segundo pedaço
#########
ggmonthplot(selic_b) +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a - média)") +
labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")
# Gráfico por ano e mês
ggseasonplot(selic_b, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")
m <- seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
plot(m)
summary(m)
##
## Call:
## seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AR-Nonseasonal-01 -0.19030 0.14589 -1.304 0.192
## AR-Nonseasonal-02 0.09671 0.09194 1.052 0.293
## AR-Nonseasonal-03 0.51926 0.08351 6.218 5.03e-10 ***
## MA-Nonseasonal-01 0.12082 0.17475 0.691 0.489
## MA-Seasonal-12 0.99913 0.12930 7.727 1.10e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (3 1 1)(0 1 1) Obs.: 87 Transform: none
## AICc: 306.6, BIC: 319.2 QS (no seasonality in final):3.007
## Box-Ljung (no autocorr.): 44.68 ** Shapiro (normality): 0.9865
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
## been found in one or more of the estimated spectra.
qs(m)
## qs p-val
## qsori 30.28836 0.00000
## qsorievadj 30.28836 0.00000
## qsrsd 0.00000 1.00000
## qssadj 3.00726 0.22232
## qssadjevadj 3.00726 0.22232
## qsirr 0.00000 1.00000
## qsirrevadj 0.00000 1.00000
summary(m)
##
## Call:
## seas(x = selic_b, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AR-Nonseasonal-01 -0.19030 0.14589 -1.304 0.192
## AR-Nonseasonal-02 0.09671 0.09194 1.052 0.293
## AR-Nonseasonal-03 0.51926 0.08351 6.218 5.03e-10 ***
## MA-Nonseasonal-01 0.12082 0.17475 0.691 0.489
## MA-Seasonal-12 0.99913 0.12930 7.727 1.10e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (3 1 1)(0 1 1) Obs.: 87 Transform: none
## AICc: 306.6, BIC: 319.2 QS (no seasonality in final):3.007
## Box-Ljung (no autocorr.): 44.68 ** Shapiro (normality): 0.9865
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
## been found in one or more of the estimated spectra.
selic_2b <- final(m)
ggmonthplot(selic_2b) +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média diária e mensal da Taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")
# Gráfico por ano e mês
ggseasonplot(selic_2b, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média mensal da Taxa Selic", x = "Mês", subtitle = "Mai/99 a Jul/06")
selic_2b <- final(m)
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1999-05-01"), length = length(selic_b), by = "1 month"), selic_b, selic_2b)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = selic_b, colour = "Original")) +
geom_line(aes(y = selic_2b, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Selic (%a.a.)")+
labs(color="Série") +
theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1
#############
# Terceiro pedaço
#############
ggmonthplot(selic_c) +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a - média)") +
labs(title = "Média diária e mensal da taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")
# Gráfico por ano e mês
ggseasonplot(selic_c, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média mensal da taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")
m <- seas(x = selic_c, transform.function = "none", regression.aictest = NULL)
plot(m)
summary(m)
##
## Call:
## seas(x = selic_c, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AR-Nonseasonal-01 -0.15177 0.10739 -1.413 0.157587
## AR-Nonseasonal-02 0.23574 0.07095 3.323 0.000891 ***
## AR-Nonseasonal-03 0.58420 0.06552 8.916 < 2e-16 ***
## MA-Nonseasonal-01 0.31636 0.12889 2.455 0.014107 *
## MA-Seasonal-12 0.99998 0.08835 11.319 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (3 1 1)(0 1 1) Obs.: 136 Transform: none
## AICc: 305.7, BIC: 321.9 QS (no seasonality in final): 2.96
## Box-Ljung (no autocorr.): 57.69 *** Shapiro (normality): 0.985
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
## been found in one or more of the estimated spectra.
qs(m)
## qs p-val
## qsori 44.66623 0.00000
## qsorievadj 44.66623 0.00000
## qsrsd 0.00000 1.00000
## qssadj 2.96000 0.22764
## qssadjevadj 2.96000 0.22764
## qsirr 0.00000 1.00000
## qsirrevadj 0.00000 1.00000
## qssori 33.30455 0.00000
## qssorievadj 33.30455 0.00000
## qssrsd 0.00000 1.00000
## qsssadj 3.39334 0.18329
## qsssadjevadj 3.39334 0.18329
## qssirr 0.00000 1.00000
## qssirrevadj 0.00000 1.00000
selic_2c <- final(m)
ggmonthplot(selic_2c) +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média diária e mensal da Taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")
# Gráfico por ano e mês
ggseasonplot(selic_2c, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Selic (%a.a. - média)") +
labs(title = "Média mensal da Taxa Selic", x = "Mês", subtitle = "Ago/06 a Nov/17")
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("2006-08-01"), length = length(selic_c), by = "1 month"), selic_c, selic_2c)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = selic_c, colour = "Original")) +
geom_line(aes(y = selic_2c, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Selic (%a.a.)")+
labs(color="Série") +
theme_bw() + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 9))
p1
###############
# Juntando novamente as três séries
###############
comb <- ts.union(selic_2a, selic_2b, selic_2c)
selic_final <- pmin(comb[,1], comb[,2], comb[,3], na.rm = TRUE)
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(selic_final), by = "1 month"), selic, selic_final)
names(df) <- c("Data", "Original", "Série Filtrada")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = selic, colour = "Original")) +
geom_line(aes(y = selic_final, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Selic (%a.a.)")+
labs(color="Série") +
theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1
m <- seas(x = selic, transform.function = "none", regression.aictest = NULL)
summary(m)
##
## Call:
## seas(x = selic, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## LS1996.Mar -4.09295 0.97485 -4.199 2.69e-05 ***
## LS1997.Nov 21.72884 0.95722 22.700 < 2e-16 ***
## LS1998.Jan -6.54038 0.99031 -6.604 3.99e-11 ***
## LS1998.Feb -5.49776 0.98955 -5.556 2.76e-08 ***
## LS1998.Apr -7.40308 0.95723 -7.734 1.04e-14 ***
## LS1998.Sep 15.32547 1.03446 14.815 < 2e-16 ***
## AO1998.Oct 7.50672 0.97560 7.694 1.42e-14 ***
## AO1998.Nov 3.92048 0.85822 4.568 4.92e-06 ***
## AO1999.Jan -4.09824 0.92123 -4.449 8.64e-06 ***
## AO1999.Mar 14.93677 0.85625 17.444 < 2e-16 ***
## LS1999.May -5.36681 1.05880 -5.069 4.00e-07 ***
## LS1999.Jun -4.48989 0.92662 -4.845 1.26e-06 ***
## AR-Nonseasonal-01 -0.25889 0.10018 -2.584 0.00976 **
## AR-Nonseasonal-02 0.13521 0.06290 2.150 0.03159 *
## AR-Nonseasonal-03 0.49771 0.05109 9.742 < 2e-16 ***
## MA-Nonseasonal-01 0.14349 0.11460 1.252 0.21053
## MA-Seasonal-12 1.00000 0.04180 23.923 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (3 1 1)(0 1 1) Obs.: 263 Transform: none
## AICc: 799.2, BIC: 859.7 QS (no seasonality in final):2.135
## Box-Ljung (no autocorr.): 59.74 *** Shapiro (normality): 0.9856 **
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
## been found in one or more of the estimated spectra.
##
## Notes:
## - Unable to test LS1998.Jan due to regression matrix
## singularity.
## - Unable to test LS1999.Jun due to regression matrix
## singularity.
## - Unable to test LS1998.Feb due to regression matrix
## singularity.
## - Unable to test LS1998.Feb due to regression matrix
## singularity.
## - Unable to test AO1998.Jan due to regression matrix
## singularity.
## - Unable to test LS1999.May due to regression matrix
## singularity.
selic2 <- final(m)
df <- data.frame(seq(as.Date("1996-01-01"), length = length(capital_trabalho2), by = "1 month"), selic, selic2, selic_final)
names(df) <- c("Data", "selic", "selic2", "selic_final")
p1 <- ggplot(df, aes(Data)) +
geom_line(aes(y = selic, colour = "Original")) +
geom_line(aes(y = selic2, colour = "Automático")) +
geom_line(aes(y = selic_final, colour = "Manual")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "Selic (%a.a.)")+
labs(color="Método") +
theme_bw() + theme(axis.text.x = element_text(angle=30, hjust = 1, size = 7))
p1
#######################
# IBC-Br #
#######################
m <- seas(x = ibcbr)
#plot(m)
#summary(m)
# Usa a interface gráfica
#view(m)
# Como ficou igual, também não mexi nele.
#######################
# PIB Mensal #
#######################
cores <- brewer.pal(6, "Dark2")
ggmonthplot(PIB) +
theme_bw() +
scale_y_continuous(name = "PIB (log-dif - média)") +
labs(title = "Média diária e mensal da variação do PIB", x = "Mês", subtitle = "Nenhuma transformação")
# Gráfico por ano e mês
ggseasonplot(PIB, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "PIB (log-dif - média)") +
labs(title = "Média mensal da da variação do PIB", x = "Mês", subtitle = "Nenhuma transformação")
m <- seas(x = PIB, transform.function = "none", regression.aictest = NULL)
## Model used in SEATS is different: (2 0 2)(0 1 1)
plot(m)
summary(m)
##
## Call:
## seas(x = PIB, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AR-Nonseasonal-01 -1.10213 0.06988 -15.771 < 2e-16 ***
## AR-Nonseasonal-02 -0.38179 0.08846 -4.316 1.59e-05 ***
## AR-Nonseasonal-03 -0.13937 0.06474 -2.153 0.0313 *
## MA-Nonseasonal-01 -0.96137 0.03241 -29.666 < 2e-16 ***
## MA-Seasonal-12 0.79111 0.04043 19.568 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (3 0 1)(0 1 1) Obs.: 263 Transform: none
## AICc: -1192, BIC: -1171 QS (no seasonality in final):0.8726
## Box-Ljung (no autocorr.): 86.12 *** Shapiro (normality): 0.9961
## Messages generated by X-13:
## Warnings:
## - At least one visually significant trading day peak has
## been found in one or more of the estimated spectra.
##
## Notes:
## - Model used for SEATS decomposition is different from the
## model estimated in the regARIMA modeling module of
## X-13ARIMA-SEATS.
# Usa a interface gráfica
#view(m)
# Salva a série final em uma nova variável
PIB2 <- final(m)
# Compara a série ajustada com a série original
df <- data.frame(seq(as.Date("1996-01-01"), length = length(PIB), by = "1 month"), PIB, PIB2)
names(df) <- c("Data", "Original", "Série Filtrada")
p2 <- ggplot(df, aes(Data)) +
geom_line(aes(y = PIB, colour = "Original")) +
geom_line(aes(y = PIB2, colour = "Filtrada")) +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%m-%Y"))+
scale_y_continuous(name = "PIB Mensal (variação mensal do índice, base = 11/17)")+
labs(color="Variável") +
theme_bw()
p2 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 8))
#######################
# Câmbio #
#######################
cores <- brewer.pal(6, "Dark2")
ggmonthplot(cambio) +
theme_bw() +
scale_y_continuous(name = "Tx. Cambio (log-dif - média)") +
labs(title = "Média diária e mensal da variação do Câmbio", x = "Mês", subtitle = "Nenhuma transformação")
# Gráfico por ano e mês
ggseasonplot(PIB, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "Tx. Cambio (log-dif - média)") +
labs(title = "Média mensal da variação do Câmbio", x = "Mês", subtitle = "Nenhuma transformação")
m <- seas(x = cambio, transform.function = "none", regression.aictest = NULL)
plot(m)
summary(m) # Summary diz que não deve fazer ajuste sazonal
##
## Call:
## seas(x = cambio, transform.function = "none", regression.aictest = NULL)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## AO1999.Jan 0.48976 0.04001 12.241 < 2e-16 ***
## AO1999.Mar -0.18202 0.04001 -4.549 5.38e-06 ***
## AO2002.Jul 0.18745 0.04001 4.685 2.80e-06 ***
## AO2002.Sep 0.27601 0.04001 6.899 5.25e-12 ***
## AO2011.Sep 0.16420 0.04001 4.104 4.06e-05 ***
## AR-Nonseasonal-01 0.11781 0.06123 1.924 0.0544 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## SEATS adj. ARIMA: (1 0 0) Obs.: 263 Transform: none
## AICc: -928.6, BIC: -904 QS (no seasonality in final):0.4173
## Box-Ljung (no autocorr.): 31.84 Shapiro (normality): 0.9713 ***
## Messages generated by X-13:
## Warnings:
## - Series should not be a candidate for seasonal adjustment
## because the spectrum of the prior adjusted series (Table
## B1) has no visually significant seasonal peaks.
# Usa a interface gráfica
#view(m)
#######################
# IPCA #
#######################
cores <- brewer.pal(6, "Dark2")
ggmonthplot(ipca) +
theme_bw() +
scale_y_continuous(name = "IPCA (% Acum. 12 meses - média)") +
labs(title = "Média diária e mensal do IPCA acum. 12 meses", x = "Mês", subtitle = "Nenhuma transformação")
# Gráfico por ano e mês
ggseasonplot(ipca, year.labels = TRUE) +
geom_point() +
theme_bw() +
scale_y_continuous(name = "IPCA (% Acum. 12 meses - média)") +
labs(title = "Média mensal do IPCA acum. 12 meses", x = "Mês", subtitle = "Nenhuma transformação")
m <- seas(x = ipca, transform.function = "none", regression.aictest = NULL)
## Model used in SEATS is different: (0 2 1)
#plot(m)
#summary(m) # Series should not be a candidate for seasonal adjustment because the spectrum of the prior adjusted series (Table B1) has no visually significant seasonal peaks.
# Usa a interface gráfica
#view(m)
# Repete o mesmo gráfico de antes, mas com as variáveis dessazonalizadas
df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho2, selic2, ibcbr, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "IBCBr", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")
# Gráficos individuais
p1 <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Capital trabalho") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p1 <- p1 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p2 <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[2])+
scale_y_continuous(name="Selic (%a.a.)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p2 <- p2 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p3 <- ggplot(df2[which(df2$variable == "IBCBr"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[3])+
scale_y_continuous(name="IBC-Br") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p3 <- p3 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p4 <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[4])+
scale_y_continuous(name="Tx. Câmbio (var)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p4 <- p4 + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p5 <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[5])+
scale_y_continuous(name="IPCA (acum. 12m.)") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
theme_bw()
p5 <- p5 + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
grid.arrange(p1, p2, p3, p4, p5, ncol=1, nrow = 5)
# Repete o mesmo gráfico de antes, mas com as variáveis dessazonalizadas
cores <- brewer.pal(6, "Dark2")
df1 <- data.frame(seq(as.Date(inicio), length = length(ipca), by = "1 month"), capital_trabalho_final, selic_final, PIB2, cambio, ipca)
names(df1) <- c("Data", "Capital_trabalho", "Selic", "PIB", "Cambio", "IPCA")
df2 <- melt(data = df1, id.vars = "Data")
## Warning: attributes are not identical across measure variables; they will
## be dropped
# Gráficos individuais
p1a <- ggplot(df2[which(df2$variable == "Capital_trabalho"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[1])+
scale_y_continuous(name="Capital trabalho") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p1a <- p1a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p2a <- ggplot(df2[which(df2$variable == "Selic"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[2])+
scale_y_continuous(name="Selic (%a.a.)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p2a <- p2a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p3a <- ggplot(df2[which(df2$variable == "PIB"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[3])+
scale_y_continuous(name="PIB") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p3a <- p3a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p4a <- ggplot(df2[which(df2$variable == "Cambio"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[4])+
scale_y_continuous(name="Tx. Câmbio (var)") +
scale_x_date(date_breaks = "12 months")+
theme_bw()
p4a <- p4a + theme(axis.text.x=element_blank(), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
p5a <- ggplot(df2[which(df2$variable == "IPCA"),], aes(Data, value, colour = variable)) +
geom_line(alpha = 1, show.legend=F, colour = cores[5])+
scale_y_continuous(name="IPCA (acum. 12m.)") +
scale_x_date(date_breaks = "12 months", name = "Data", labels = date_format("%Y"))+
theme_bw()
p5a <- p5a + theme(axis.text.x = element_text(angle=25, hjust = 1, size = 6), axis.title.x = element_blank(), axis.title.y = element_text(size = 6))
grid.arrange(p1a, p2a, p3a, p4a, p5a, ncol=1, nrow = 5)
Parece que o X-13ARIMA fez um trabalho melhorzinho, então vou continuar com a série capital_trabalho_adj.
Fiz as descritivas para a série original e depois para a série trabalhada manualmente e pelo X-13ARIMA. Fiz as descritivas das variáveis que vão entrar no VAR. Depois tem que lembrar de fazer das que efetivamente ficaram Uma ideia é colocar naqueles arquivos dos VAR individuais. Acho que vou fazer isso depois…
### Descriptives
descriptives <- matrix(NA, nrow = 8, ncol = (ncol(df1)-1))
rownames(descriptives) <- c("Observações", "Mínimo", "1o quartil",
"Média", "Mediana", "3o quartil", "Máximo",
"Desv. Pad.")
colnames(descriptives) <- names(df1)[-1]
desc <- function(x) {
n <- length(x)
minimum <- min(x, na.rm = TRUE)
first_q <- quantile(x, 0.25, na.rm = TRUE)
media <- mean(x, na.rm = TRUE)
mediana <- median(x, na.rm = TRUE)
third_q <- quantile(x, 0.75, na.rm = TRUE)
maximum <- max(x, na.rm = TRUE)
std <- sd(x, na.rm = TRUE)
return(list(n = n, minimum = minimum, first_quar = first_q, media = media, mediana = mediana, third_quar = third_q, maximum = maximum, std = std))
}
for (i in 1:8){
descriptives[i, 1] <- round(as.numeric(desc(df1[,2])[i]),4)
descriptives[i, 2] <- round(as.numeric(desc(df1[,3])[i]),4)
descriptives[i, 3] <- round(as.numeric(desc(df1[,4])[i]),4)
descriptives[i, 4] <- round(as.numeric(desc(df1[,5])[i]),4)
descriptives[i, 5] <- round(as.numeric(desc(df1[,6])[i]),4)
}
descriptives[1,] <- as.integer(descriptives[1,])
descriptives <- data.frame(descriptives)
#stargazer(descriptives, summary=FALSE, header = TRUE, type = 'html')
stargazer(descriptives, summary=FALSE, header = TRUE, type = 'latex')
% Table created by stargazer v.5.2 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu % Date and time: qua, jan 24, 2018 - 03:16:31
O pacote que tem o modelo implementado é o bvarsv. Aqui (https://github.com/FK83/bvarsv/blob/master/bvarsv_Nov2015_website.pdf) e aqui (https://github.com/FK83/bvarsv/blob/master/bvarsv_replication.pdf) tem exemplos de uso.
Especificações da função bvar.sv.tvp: